Thermal conductivity of single layer MoS₂

Nisarg Vyas

School of Physical Sciences, National Institute of Science Education and Research, Jatni, Odisha 752050

Prof. Prabal Maiti

Department of Physics, Indian Institute of Science, Bengaluru 560012

Abstract

A popular representative of transition metal dichalcogenides, MoS₂ is taken up for studying the thermal properties and their dependence on temperature. The material is considered only in its 2-dimensional form with a single layer. The atomic interactions are modelled using Stillinger-Weber empirical type potential. We compute values for the thermal conductivity of mono-layer MoS₂ via Green-Kubo formalism based equilibrium molecular dynamics method. The values thus obtained, match with other computational works. In the presented work, it is also verified that the thermal conductivity of the considered material scales inversely with the system temperature. All of the molecular dynamics simulations are performed in large-scale atomic/molecular massively parallel simulator (LAMMPS).

Keywords: transition metal dichalcogenides, thermal conductivity, equilibrium molecular dynamics, Green-Kubo, Stillinger-Weber, LAMMPS

Abbreviations

INTRODUCTION

Background

Graphene, owing to its excellent thermal conductivity and high intrinsic mobility, has been a toy tool of theoretical physicists and has been studied extensively in recent years. The peculiar way in which the thermal conductivity of graphene scales as one goes from bulk to few-layer and to mono-layer graphene also makes it an interesting object to study thermal properties of matter and to dwell deeper into the mechanisms of heat conduction phenomena. However, due to the lack of electronic band gap in graphene, a very limited number of applications in nano electronics have been possible.

To overcome this deficiency, researchers have shifted their attention to another class of 2-dimensional materials: transition metal dichalcogenides (TMDCs). These materials are semiconducting in nature and can be represented as MX₂ where M can be a transition metal (like Mo, W etc.) and X can be a chalcogen (S, Se or Te). These semiconductors possess rich electronic and photonic properties with promises for novel applications in nanoelectronics. One of the applications of TMDCs is in thermo electronics where strict control over thermal properties is required. For effective implementation of the materials, it is necessary to gain a better understanding of thermal transport in them.

Objectives of the Research

In literature, the computationally predicted thermal conductivity values for single layer MoS₂ vary from 1.35-5.8 W/(m.K) [​Liu, et al, 2013​] to 103 W/(m.K) [​Gu, et al, 2014​]. The largely underestimated value obtained earlier by Liu et al, were due to the force field used. However the later studies done by Li, et al, 2013 using ab-initio methods and Boltzmann transport equations give an estimate of 83 W/(m.K). This is quite close to the experimentally obtained values of 84 ± 17 W/(m.K), achieved by Zhang, et al, 2015 using opto-thermal Raman technique. The objectives taken up in the currently undertaken study is to calculate the in-plane thermal conductivity value for mono-layer MoS₂ using SW potential and compare the results with the existing experimental values and other simulations done either with the same potential or different ones. The dependence of thermal conductivity on temperature is also studied between 50 and 450 K.

LITERATURE REVIEW

The Force Field

To properly account for the structural, vibrational and mechanical properties of the material, the Stillinger-Weber (SW) potential developed by Cem Sevik, 2016 is used. It has been shown that the interatomic interactions described by SW type empirical potentials give thermal properties that match the experimentally calculated valuues [Yan, et al, 2014]. It has also been noted that the parameters produced in this potential are in accordance with results from ab-initio methods of density functional theory for both MoS₂ and MoSe₂ [Gu, et al, 2014].

The SW potential as proposed in 1985 [Stillinger, et al, 1985] defines on site potential as:

​$U_i=\frac{1}{2}\sum _{i\ne j}{V}_2\left(r_{ij}\right)+\frac{1}{2}\sum _{j\ne i}\sum _{k\ne i,j}{h}_{ijk}$

Where V₂ describes the two body terms called the bond-stretching terms, while h_ijk defines three body terms appropriately called bond-bending terms. Note that two-index terms work for a pair of atoms while three-index quantities are detrimental for triplets of atoms. The individual potential form can be given in terms of material specific parameters $A_{ij},B_{ij},{\epsilon }_{ij},{\sigma }_{ij},a_{ij},{\lambda }_{ijk},{\gamma }_{ij}$and $\cos \left({\theta }_{ijk}\right)$ as follows.

$V_2\left(r_{ij}\right)=A_{ij}{\epsilon }_{ij}\left[B_{ij}{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^4-1\right]exp\left(\frac{1}{{\displaystyle \raisebox{1ex}{$r_{ij}$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{ij}$}\right.}-a_{ij}}\right)$

$h_{ijk}={\epsilon }_{ij}{\lambda }_{ijk}exp\left[\frac{{\gamma }_{ij}}{\frac{r_{ij}}{{\sigma }_{ij}}-a_{ij}}+\frac{{\gamma }_{ik}}{\frac{r_{ik}}{{\sigma }_{ik}}-a_{ik}}\right]\times {\left(\cos \left({\theta }_{ijk}\right)-\cos \left({\theta }_{0ijk}\right)\right)}^2$​

METHODOLOGY

Concepts

There are quite a few well proven methods in molecular dynamics to calculate thermal conductivity. Among these are the equilibrium molecular dynamics (EMD) method based on Green-Kubo relation and Fourier law based other non-equilibrium molecular dynamics (NEMD) methods. When the heat transport in diffusive, both of the methods give consistent calculations on thermal properties. However we use the Green-Kubo formalism in the present work to estimate the thermal conductivity of mono-layer MoS_2.

The Green-Kubo formalism works as follows. The ensemble average of heat current auto-correlation function is calculated and then is numerically integrated to obtain a value of thermal conductivity. The formula which relates the ensemble averaged quantity with, thermal conductivity, kappa is given as follows:

​$\begin{array}{l}\kappa \hspace{0.25em} =\hspace{0.25em} \frac{V}{k_B{T}^2}{\int }_0^{\mathrm{\infty }}⟨J_x(0)J_x(t)⟩dt\\ \\ \end{array}$ $\begin{array}{l}\hspace{0.25em} =\hspace{0.25em} \frac{V}{3k_B{T}^2}{\int }_0^{\mathrm{\infty }}⟨J(0)\mathrm{.}J(t)⟩dt\\ \\ \end{array}$

where k_B is the Boltzmann's constant, T is the system temperature and V is the system volume. Note that the above two expressions equate only when isotropic bulk systems are considered. In the present case of 2-D MoS₂ sheet, it being anisotropic, only the x component of heat flux is considered in calculations to give a thermal conductivity value along the length (in the direction of the x axis). The heat current J is calculated as

$J=\frac{1}{V}\left[\sum _i{e}_i{v}_i-\sum _i{S}_i{v}_i\right]$

In the first term, e_i are the energy terms (both potential and kinetic) while V_i is a 3×3 matrix multiplied to energy vector to get a vector again. The S_i is per atom stress tensor. Both of these quantites are computes at every timestep by appropriate compute commands in LAMMPS.

Methods

The system used in the study is a 11.2 nm long and 2.2 nm wide ribbon of mono-layer MoS₂. The out of plane direction is the direction of the z axis while the length of the ribbon is along the x axis (see fig. 1 and 2). The simulation box has periodic boundary conditions in the direction of the y and z axes while the box’s boundaries are non-periodic in nature along the length. However, the effect of periodic boundaries along the z axis is nullified by keeping sufficient vacuum on both sides of the ribbon. The initial structure prepared with approximate interatomic distances is minimized for the potential energy using the conjugate gradient algorithm implemented in LAMMPS. After the minimization step, the correctness of the MoS₂ structure is determined by observing the radial distribution function with OVITO, a visualization software. Peaks seen in the radial distribution function, which shows the variation of density as a function of distance from a particle, appear at values that correspond to the nearest and second nearest neighbour distance for MoS₂.

Fig 1 Out of plane or top view of the MoS₂ ribbon. White spheres are Sulphur atoms while orange ones are Molybdenum atom.

a) Edge view along the y axis (along the width)

b) Side view along the x axis (along the length). Note that the sheet ends in two different ways with Mo on one side and with S atoms on the other.

Fig 2 Structure of MoS₂ sheet used in the simulation. Visualization in OVITO

The protocol followed throughout the simulations is as follows. The atoms are given a Maxwellian velocity distribution in accordance with the desired temperature for the system. The prepared structure is then kept in a canonical ensemble (nvt) for 500 ps to equilibrate the temperature of the system. Once the system reaches a steady state, the heat current at every 10^th step is used for calculation for 5×10⁶ steps. This gives a numerical integration of the heat current auto-correlation function over a total time period of 5 ns. While the integral is evaluated, the system simultaneously remains in a micro-canonical ensemble (nve). Note that the Green-Kubo method is called EMD because the heat current is calculated when the system resides in a steady equilibrium state, this justifies the invoking of micro-canonical ensemble.

Fig 3 On x axis is the number of simulation steps versus thermal conductivity values to show convergence.

The kappa calculation and averaging is done only after making sure that the auto-correlated flux values have converged. See fig. 3 for similar check made on the calculation for system at 300 K. The volume of the system, needed to calculate the thermal conductivity using equation in ​Sec. 3.1​, is computed using surface mesh analysis in OVITO which computes the solid volume after constructing a closed geometrical surface enclosing all the atoms within it [Stukowski, et al, 2013]. For simple systems like the one considered here, the method is doubtlessly accurate.

RESULTS AND DISCUSSION

After equilibrating the system at different temperatures, each time the thermal conductivity calculations are averaged for the last 1 ns of the simulation time. The average value when plotted against temperature with the standard deviation of the data as the error bar, gives the following plot (fig. 4).

Fig 4 The fit obtained for k vs. T; f(x) = 84.25*(1/T) + 6.89

The thermal conducivity estimation, for a MoS₂ sheet of approximate length 11 nm, at 300 K is around 7.1 W/(m.K) which lies close to 5 W/(m.K), an estimation given by Hong, et al, 2016 at the same length and temperature and using the same force fields. However their method of estimation is Green-Kubo as well as non-equilibrium molecular dynamics method based on Fourier's law. The 1/T curve fitting with the corresponding kappa values has a Chi square value of 0.3.

CONCLUSION

The obtained values of thermal conductivity for single layer MoS₂ are close to those obtained by Hong, et al, 2016. However, it will be too hasty to deduce the inverse dependence of kappa on the system temperature as the plot does not have satisfactory number of data points. We realize that the claim that kappa depends inversely on temperature can be strengthened only by generating statistically enough data points. Increasing the number of simulation runs per data point would not only give more reliable thermal conductivity values but would also prevent the underestimation of errors.

Nevertheless, the effort to compute and cross check the thermal conductivity of single layer MoS₂ can be safely said to have been accomplished. Further studies in this direction can be taken up with confidence. Similar to the heterostructure proposed by Yang, et al, 2012, it has not escaped our notice that a heterostructure of MoS₂ and MoSe₂ can be exploited to achieve a material of thermal conductivity lower than both of the materials individually. This can be taken up as a further study steming from the one done here.

ACKNOWLEDGEMENTS

No man is an island. In this journey of about two months I have been imensely fortunate to work and interact with good people who helped me in various ways and by various means. Thanks to Vignesh for being a strict mentor. This work would not have completed without you. I would like to extend heartiest thanks to Shubhadeep Dasgupta and Jeet Majumdar for their constant support and for the insights they provided into methods of molecular dynamics, without which this work would not have made much sense to me. A sincere thanks goes to Indrajit Maity for helping me with his expertize whenever I was in trouble. A cheerful thanks to my cousins, Jayaswini and Sumedha, and to my uncle and aunt for making my stay in Bangalore so comfortable and for making my weekends enjoyful!

Lastly, I wish to express my gratitude towards Prof. Prabal K Maiti and the Indian Academy of Sciences for providing me this opportunity to learn and to grow.